function [betaXhr,betaXhrRegime,R2_rec,R2_exp,nameRegressors] = forecastRegression(yields,shortRate,maturities,REC,h)
% Data inputs:
% yields: T * numYields where yields are annually and continuously compounded
% maturities: 1*numYields with maturities in months
% RECdummy: 0 1 regression dummy, where 0 is boom and 1 is recression
% h: sampletime_per_maturityTime

% Preable
T         = size(yields,1);
numYields = length(maturities);

% forwards: f[t]_(n-h)_n = log(P[t]_(n-h)_n) - log(P[t]_n_n)
%                     = -(n-h)*y[t]_(n-h) - (-(n)*y[t]_(n))
%                     = n*y[t]_(n)        -(n-h)*y[t]_(n-h)
% where y(n) is the spot yield at maturity n
% h-month forward spreads
annualMaturities = maturities/12;
forward  = repmat(annualMaturities(1,2:end),T,1)  .*yields(:,2:end) ...
          -repmat(annualMaturities(1,1:end-1),T,1).*yields(:,1:end-1);
if T < 1000
    % I.e. we have empirical data and we smooth the forward rates
    smoother.bw=6;
    smoother.kernel=1;
    forward = smoothingFct(forward,maturities(1,2:end)', smoother );
end
    
% The forward spread: forward rate - baseline yield.
% Note that the baseline yield is in annual terms so we divide by 12 and
% multiply by h
%forwardSpreads   = forward-repmat(yields(:,1),1,numYields-1)*h/12;          
forwardSpreads   = forward-repmat(shortRate,1,numYields-1)*h/12;          

% Term spread: y[t](n)-y[t](1) for all n.
%termSpreads     = yields(:,2:end)-repmat(yields(:,1),1,numYields-1);
termSpreads     = yields(:,2:end)-repmat(shortRate,1,numYields-1);

% Ex post holding period returns: hr[t]_n = n*y[t-h](n) - (n-h)*y[t](n-h)
holdPerReturn   = repmat(annualMaturities(1,2:end),T-h,1)  .*yields(1:end-h,2:end) ...
                   -repmat(annualMaturities(1,1:end-1),T-h,1).*yields(1+h:end,1:end-1);

% Ex post excess holding period returns: xhr[t]_n = hr[t]_t - y[t-h](h)*h/12
%excessHoldPerReturn = holdPerReturn - repmat(yields(1:end-h,1),1,numYields-1)*h/12;               
excessHoldPerReturn = holdPerReturn - repmat(yields(1:end-h,maturities == h),1,numYields-1)*h/12;               
               
%% Generating the CP factor, computed based on one-year forward rates
[~,CP_idx]=find(repmat(annualMaturities,5,1)==repmat([1:5]',1,numYields));
forward1Year  = repmat(annualMaturities(1,CP_idx(2:end)),T,1)  .*yields(:,CP_idx(2:end)) ...
               -repmat(annualMaturities(1,CP_idx(1:end-1)),T,1).*yields(:,CP_idx(1:end-1));
holdPerReturn1Year = repmat(annualMaturities(1,CP_idx(2:end)),T-12,1)  .*yields(1:end-12,CP_idx(2:end)) ...
                    -repmat(annualMaturities(1,CP_idx(1:end-1)),T-12,1).*yields(1+12:end,CP_idx(1:end-1));
excessholdPerReturn1Year  = holdPerReturn1Year - repmat(yields(1:end-12,CP_idx(1)),1,4);
               
Y     = mean(excessholdPerReturn1Year,2);
X     = [ones(T-12,1) yields(1:end-12,CP_idx(1)) forward1Year(1:end-12,:)];
gamma = olsgmm(Y,X,12,1);
CP    = [ones(T,1)    yields(1:end,CP_idx(1))    forward1Year(1:end,:)]*gamma;

%% Runing the forecast regressions for all maturities and for the three predictors
index_RECh    = find(REC(1:end-h,1)==1);
index_EXPh    = find((1-REC(1:end-h,1))==1);
betaXhr       = NaN(3,numYields);
betaXhrRegime = NaN(3,2,numYields);
R2_rec        = NaN(3,numYields);
R2_exp        = NaN(3,numYields);
startMat      = find(maturities >= 2*12);
nameRegressors= {'Term spread','Forward spread','CP factor'};
for j=1:3
    %for n=1:numYields-1
    for n=startMat:numYields-1
        % excessHoldPerReturn regressed on term spread
        Y = excessHoldPerReturn(:,n);
        if j == 1
            X = [ones(size(Y,1),1) termSpreads(1:end-h,n)];
        elseif j == 2
            X = [ones(size(Y,1),1) forwardSpreads(1:end-h,n)];
        elseif j == 3
            X = [ones(size(Y,1),1) CP(1:end-h,1)];
        end
        
        tmp = olsgmm(Y,X,h,1);
        betaXhr(j,n)   = tmp(2,1);
        residuals      = Y - X*tmp;
        residualsConst = Y - mean(Y);
        
        % Regime dependent: excessHoldPerReturn regressed on term spread
        X1=[X.*repmat(1-REC(1:end-h,1),1,2) X.*repmat(REC(1:end-h,1),1,2) ];
        tmp = olsgmm(Y,X1,h,1);
        betaXhrRegime(j,:,n) = [tmp(2,1);tmp(4,1)];
        
        R2_rec(j,n) = (1-sum(residuals(index_RECh).^2)/sum(residualsConst(index_RECh).^2))*100;
        R2_exp(j,n) = (1-sum(residuals(index_EXPh).^2)/sum(residualsConst(index_EXPh).^2))*100;
    end
end

end

